/*!
 * \file
 * \brief Implementation of modified Bessel functions of order one
 * \author Tony Ottosson
 *
 * -------------------------------------------------------------------------
 *
 * Copyright (C) 1995-2010  (see AUTHORS file for a list of contributors)
 *
 * This file is part of IT++ - a C++ library of mathematical, signal
 * processing, speech processing, and communications classes and functions.
 *
 * IT++ is free software: you can redistribute it and/or modify it under the
 * terms of the GNU General Public License as published by the Free Software
 * Foundation, either version 3 of the License, or (at your option) any
 * later version.
 *
 * IT++ is distributed in the hope that it will be useful, but WITHOUT ANY
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 * details.
 *
 * You should have received a copy of the GNU General Public License along
 * with IT++.  If not, see <http://www.gnu.org/licenses/>.
 *
 * -------------------------------------------------------------------------
 *
 * This is slightly modified routine from the Cephes library:
 * http://www.netlib.org/cephes/
 */

#include <itpp/base/bessel/bessel_internal.h>


/*
 * Modfied Bessel function of order one
 *
 * double x, y, i1();
 *
 * y = i1( x );
 *
 * DESCRIPTION:
 *
 * Returns modified Bessel function of order one of the
 * argument.
 *
 * The function is defined as i1(x) = -i j1( ix ).
 *
 * The range is partitioned into the two intervals [0,8] and
 * (8, infinity).  Chebyshev polynomial expansions are employed
 * in each interval.
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0, 30       30000       1.9e-15     2.1e-16
 */

/*
 * Modified Bessel function of order one, exponentially scaled
 *
 * double x, y, i1e();
 *
 * y = i1e( x );
 *
 * DESCRIPTION:
 *
 * Returns exponentially scaled modified Bessel function
 * of order one of the argument.
 *
 * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0, 30       30000       2.0e-15     2.0e-16
 * See i1().
 */

/*
  Cephes Math Library Release 2.8:  June, 2000
  Copyright 1985, 1987, 2000 by Stephen L. Moshier
*/


/* Chebyshev coefficients for exp(-x) I1(x) / x
 * in the interval [0,8].
 *
 * lim(x->0){ exp(-x) I1(x) / x } = 1/2.
 */

static double A[] = {
  2.77791411276104639959E-18,
  -2.11142121435816608115E-17,
  1.55363195773620046921E-16,
  -1.10559694773538630805E-15,
  7.60068429473540693410E-15,
  -5.04218550472791168711E-14,
  3.22379336594557470981E-13,
  -1.98397439776494371520E-12,
  1.17361862988909016308E-11,
  -6.66348972350202774223E-11,
  3.62559028155211703701E-10,
  -1.88724975172282928790E-9,
  9.38153738649577178388E-9,
  -4.44505912879632808065E-8,
  2.00329475355213526229E-7,
  -8.56872026469545474066E-7,
  3.47025130813767847674E-6,
  -1.32731636560394358279E-5,
  4.78156510755005422638E-5,
  -1.61760815825896745588E-4,
  5.12285956168575772895E-4,
  -1.51357245063125314899E-3,
  4.15642294431288815669E-3,
  -1.05640848946261981558E-2,
  2.47264490306265168283E-2,
  -5.29459812080949914269E-2,
  1.02643658689847095384E-1,
  -1.76416518357834055153E-1,
  2.52587186443633654823E-1
};


/* Chebyshev coefficients for exp(-x) sqrt(x) I1(x)
 * in the inverted interval [8,infinity].
 *
 * lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi).
 */

static double B[] = {
  7.51729631084210481353E-18,
  4.41434832307170791151E-18,
  -4.65030536848935832153E-17,
  -3.20952592199342395980E-17,
  2.96262899764595013876E-16,
  3.30820231092092828324E-16,
  -1.88035477551078244854E-15,
  -3.81440307243700780478E-15,
  1.04202769841288027642E-14,
  4.27244001671195135429E-14,
  -2.10154184277266431302E-14,
  -4.08355111109219731823E-13,
  -7.19855177624590851209E-13,
  2.03562854414708950722E-12,
  1.41258074366137813316E-11,
  3.25260358301548823856E-11,
  -1.89749581235054123450E-11,
  -5.58974346219658380687E-10,
  -3.83538038596423702205E-9,
  -2.63146884688951950684E-8,
  -2.51223623787020892529E-7,
  -3.88256480887769039346E-6,
  -1.10588938762623716291E-4,
  -9.76109749136146840777E-3,
  7.78576235018280120474E-1
};


double i1(double x)
{
  double y, z;

  z = fabs(x);
  if (z <= 8.0) {
    y = (z / 2.0) - 2.0;
    z = chbevl(y, A, 29) * z * exp(z);
  }
  else {
    z = exp(z) * chbevl(32.0 / z - 2.0, B, 25) / sqrt(z);
  }
  if (x < 0.0)
    z = -z;
  return(z);
}

double i1e(double x)
{
  double y, z;

  z = fabs(x);
  if (z <= 8.0) {
    y = (z / 2.0) - 2.0;
    z = chbevl(y, A, 29) * z;
  }
  else {
    z = chbevl(32.0 / z - 2.0, B, 25) / sqrt(z);
  }
  if (x < 0.0)
    z = -z;
  return(z);
}
